## ---- spatial-analysis ----

# Lev 1 (comune level)

library(sp)
library(rgdal)
library(sf)
library(rgeos)
library(DBI)
library(RPostgreSQL)
library(rpostgis)
library(reshape2)
library(dplyr)
library(raster)

options(scipen=999)

# Connect to local database
conn <- dbConnect(
  drv = PostgreSQL(),
  dbname = "istat_sez2011",
  host = "localhost",
  port = "5432",
  user = "francesco",
  password = "")

comuni_sp <- 
  readOGR('istat_comuni_2018_modified/', layer = 'istat_comuni_2018')

comuni2018_candidate_votes <- 
  dbReadTable(conn, c("ita_ge18", "comuni2018_camera_candidate_votes"))
comuni2018_party_votes <- 
  dbReadTable(conn, c("ita_ge18", "comuni2018_camera_party_votes"))
comuni2018_place_stats <- 
  dbReadTable(conn, c("ita_ge18", "comuni2018_camera_place_stats"))

comuni_to_merge <- c("001130", "001265", "001063", "001028", "001316")

comuni2018_candidate_votes$istat_cod[
  comuni2018_candidate_votes$istat_cod %in% comuni_to_merge] <- 
  "000000"
comuni2018_party_votes$istat_cod[
  comuni2018_party_votes$istat_cod %in% comuni_to_merge] <- 
  "000000"
comuni2018_place_stats$istat_cod[
  comuni2018_place_stats$istat_cod %in% comuni_to_merge] <- 
  "000000"

comuni2018_candidate_votes <- 
  comuni2018_candidate_votes[,c("istat_cod","coalition","candidate_votes")] %>%
  dplyr::group_by(istat_cod, coalition) %>%
  dplyr::summarise (votes = sum(candidate_votes)) %>%
  dplyr::mutate(perc = votes / sum(votes))

require(dplyr)
comuni2018_party_votes <- 
  comuni2018_party_votes[,c("istat_cod","party","party_votes")] %>%
  dplyr::group_by(istat_cod, party) %>%
  dplyr::summarise (votes = sum(party_votes)) %>%
  dplyr::mutate(perc = votes / sum(votes))

comuni2018_place_stats <- 
  comuni2018_place_stats[,c('istat_cod', "registered_voters", "effective_voters")] %>%
  dplyr::group_by(istat_cod) %>%
  dplyr::summarise(registered_voters = sum(registered_voters),
                   effective_voters = sum(effective_voters))

comuni2018_camera2013 <- dbReadTable(conn, c("ita_ge18", "comuni2018_camera2013"))

comuni_2013 <- 
  dbGetQuery(conn, "SELECT \"PRO_COM\", \"IL POPOLO DELLA LIBERTA'\",
\"LEGA NORD\",\"MOVIMENTO 5 STELLE BEPPEGRILLO.IT\",\"PARTITO DEMOCRATICO\", 
\"FRATELLI D'ITALIA\", \"ppdt_nRegisteredVoters\", 
\"ppdt_nAssignedVotes\", \"ppdt_nEffectiveVoters\" FROM elections.comuni_2013_camera 
             WHERE \"PRO_COM\" = '1130' OR \"PRO_COM\" = '1265' OR \"PRO_COM\" = '1063' 
             OR \"PRO_COM\" = '1028' OR \"PRO_COM\" = '1316';")
comuni_2013$PRO_COM <- NULL
names(comuni_2013)[1:5] <- c('PDL','LN','M5S','PD',"FdI")
comuni_2013$PRO_COM_T_2018 <- "000000"

comuni_2013 <-
  comuni_2013 %>%
  dplyr::group_by(PRO_COM_T_2018) %>%
  dplyr::summarise_all(sum)
comuni2018_camera2013 <- 
  rbind(comuni2018_camera2013, comuni_2013)

comuni_sp <- 
  merge(comuni_sp, 
        comuni2018_camera2013[,
                              c('PRO_COM_T_2018','PDL','LN','M5S','PD',
                                'ppdt_nRegisteredVoters', 'ppdt_nEffectiveVoters', 
                                'ppdt_nAssignedVotes')], 
        by.x='PRO_COM_T', by.y='PRO_COM_T_2018')
names(comuni_sp)[9:15] <- 
  paste0(names(comuni_sp)[9:15], "_2013")

comuni_sp$PDL_2013 <- comuni_sp$PDL_2013 / comuni_sp$ppdt_nAssignedVotes_2013
comuni_sp$M5S_2013 <- comuni_sp$M5S_2013 / comuni_sp$ppdt_nAssignedVotes_2013
comuni_sp$LN_2013 <- comuni_sp$LN_2013 / comuni_sp$ppdt_nAssignedVotes_2013
comuni_sp$PD_2013 <- comuni_sp$PD_2013 / comuni_sp$ppdt_nAssignedVotes_2013

comuni2018_candidate_votes <- 
  subset(comuni2018_candidate_votes, coalition == 'MOVIMENTO 5 STELLE')
comuni2018_candidate_votes <- 
  comuni2018_candidate_votes[,c('istat_cod','perc')]
names(comuni2018_candidate_votes) <- 
  c("PRO_COM_T","M5S_2018")

comuni2018_party_votes <- 
  subset(comuni2018_party_votes, party %in% c('FORZA ITALIA','LEGA','PARTITO DEMOCRATICO'))
comuni2018_party_votes <- 
  dcast(comuni2018_party_votes, istat_cod ~ party, value.var = 'perc')
names(comuni2018_party_votes) <- 
  c('PRO_COM_T','FI_2018','Lega_2018','PD_2018')

comuni2018_votes <- 
  merge(comuni2018_party_votes, comuni2018_candidate_votes, by = 'PRO_COM_T')
comuni_sp <- 
  merge(comuni_sp, comuni2018_votes, by = 'PRO_COM_T')

comuni_sp$turnout_2013 <- 
  comuni_sp$ppdt_nEffectiveVoters_2013 / comuni_sp$ppdt_nRegisteredVoters_2013

comuni_sp$ppdt_nRegisteredVoters_2013 <- NULL
comuni_sp$ppdt_nAssignedVotes_2013 <- NULL
comuni_sp$ppdt_nEffectiveVoters_2013 <- NULL

comuni2018_place_stats$turnout_2018 <- 
  comuni2018_place_stats$effective_voters / comuni2018_place_stats$registered_voters
comuni_sp <- 
  merge(comuni_sp, 
        comuni2018_place_stats[,c('istat_cod','turnout_2018')], 
        by.x = 'PRO_COM_T', by.y = 'istat_cod')

comuni_sp$M5S_diff <- (comuni_sp$M5S_2018 - comuni_sp$M5S_2013) * 100
comuni_sp$PD_diff <- (comuni_sp$PD_2018 - comuni_sp$PD_2013) * 100
comuni_sp$Lega_diff <- (comuni_sp$Lega_2018 - comuni_sp$LN_2013) * 100
comuni_sp$FI_diff <- (comuni_sp$FI_2018 - comuni_sp$PDL_2013) * 100

breaks <- c(-100, -50, -20, -10, -5, -2, 2, 5, 10, 20, 50)

for (var  in c("M5S_diff","PD_diff", "Lega_diff", "FI_diff")) {
  comuni_sp[[paste0(var, "_brk")]] <- 
    cut(comuni_sp[[var]], 
        breaks = c(breaks, max(comuni_sp[[var]], na.rm = T)),
        labels = c("-50%", "-20%" , 
                   "-10%", "-5%", 
                   "-2%", "0%", 
                   "+2%", "+5%", "+10%",
                   "+20%", "+50%"))
}

assignMacroRegion <- function(x) {
  if (x %in% c('1','3',
               '7',"2")) {
    return('North-west')
  }
  if(x %in% c('5','6',
              "4",
              "8")) {
    return('North-east')
  }
  if (x %in% c("9", "10", "11", "12")) {
    return('Centre')
  }
  if (x %in% c("13", "16", "17", 
               "15", "18", "14")) {
    return('South')
  }
  if (x %in% c("19", "20")) {
    return("Islands")
  }
}

comuni_sp$macro_region <- 
  sapply(as.character(comuni_sp$COD_REG), assignMacroRegion)
comuni_sp$macro_region <- 
  factor(comuni_sp$macro_region, 
         levels = c("North-west", "North-east",
                    "Centre", "South", "Islands" ))

res_geo_details <- 
  dbGetQuery(conn, 'SELECT "PRO_COM_T",
             ST_Area(geom) AS area
             FROM ita_ge18.comuni2018_geo;')

res_geo_details$PRO_COM_T[res_geo_details$PRO_COM_T%in% comuni_to_merge] <- "000000"
res_geo_details <- 
  res_geo_details %>%
  dplyr::group_by(PRO_COM_T) %>%
  dplyr::summarise_all(sum)

res_demographics <-
  dbGetQuery(conn, "SELECT * FROM ita_ge18.comuni2018_census2011_demographics;")

res_demographics$PRO_COM_T[
  res_demographics$PRO_COM_T %in% comuni_to_merge] <- "000000"
res_demographics$COMUNE <- NULL
res_demographics <- 
  res_demographics %>%
  dplyr::group_by(PRO_COM_T) %>%
  dplyr::summarise_all(sum)

comuni_sp <- merge(comuni_sp, res_geo_details, by = 'PRO_COM_T')
comuni_sp <- merge(comuni_sp, res_demographics, by = 'PRO_COM_T')

comuni_sp$pop_density <- comuni_sp$pop2011/comuni_sp$area
comuni_sp$unemployment <- comuni_sp$popunemployed2011/comuni_sp$popworkforce2011

comuni_sp$partwf_perc <- 
  comuni_sp$popworkforce2011 / 
  (comuni_sp$popworkforce2011 + comuni_sp$popnotworkforce2011)

comuni_sp$housewive_perc <- comuni_sp$pophousewife2011/comuni_sp$popworkforce2011
comuni_sp$over65_perc <- comuni_sp$pop65over2011/comuni_sp$pop2011
comuni_sp$degree_perc <- comuni_sp$popdegree2011/comuni_sp$pop2011

comuni_sp$foreignpop_perc <- comuni_sp$foreignpop2011/comuni_sp$pop2011
comuni_sp$foreignpop_africa_perc <- comuni_sp$foreignpopafrica2011/comuni_sp$pop2011

redditi2016 <- dbReadTable(conn, c("ita_ge18","comuni2018_redditi2016"))
comuni_sp <- merge(comuni_sp, redditi2016, by = 'PRO_COM_T', all.x = T)
comuni_sp$income2016_procapite <- comuni_sp$redditoAmmontare/comuni_sp$nContribuenti

comuni_sp$lon <- coordinates(comuni_sp)[,1]
comuni_sp$lat <- coordinates(comuni_sp)[,2]

# Lev 2 (comune level or lower)

camera_sp <- 
  readOGR('scrape_2018_ita_ge_results/output_shapefiles/', layer = 'minint_camera')

res_camera_other_candidate <- 
  dbGetQuery(conn, "SELECT * FROM ita_ge18.minint_camera_candidate_votes;")
res_camera_other_party <- 
  dbGetQuery(conn, "SELECT * FROM ita_ge18.minint_camera_party_votes;")

res_camera_other_candidate <- 
  res_camera_other_candidate %>%
  dplyr::group_by(minint_cod) %>%
  dplyr::mutate(tot_votes = sum(candidate_votes))

res_camera_other_candidate$perc <- 
  res_camera_other_candidate$candidate_votes / 
  res_camera_other_candidate$tot_votes
res_camera_other_candidate <- 
  subset(res_camera_other_candidate, 
         coalition == 'MOVIMENTO 5 STELLE',
         select = c("minint_cod", "perc"))
names(res_camera_other_candidate)[2] <- 
  'M5S_2018'

res_camera_other_party <-
  res_camera_other_party %>%
  dplyr::group_by(minint_cod) %>%
  dplyr::mutate(tot_votes = sum(party_votes))
res_camera_other_party$perc <- 
  res_camera_other_party$party_votes / res_camera_other_party$tot_votes

require(reshape2)
res_camera_other_party <- 
  dcast(data = subset(res_camera_other_party, 
                      party %in% c("LEGA", "PARTITO DEMOCRATICO", "FORZA ITALIA")), 
        formula = minint_cod ~ party, value.var = 'perc')
names(res_camera_other_party)[2:4] <- 
  c("FI_2018","Lega_2018","PD_2018")

res_place_stats <- 
  dbGetQuery(conn, "SELECT * FROM ita_ge18.minint_camera_place_stats;")
res_place_stats$turnout_2018 <- 
  res_place_stats$effective_voters / res_place_stats$registered_voters

res_camera <- merge(res_camera_other_candidate, res_camera_other_party,
                    by = 'minint_cod')
res_camera <- merge(res_camera, 
                    res_place_stats[,c('turnout_2018','minint_cod')], 
                    by = 'minint_cod')

camera_sp <- merge(camera_sp, res_camera, by = 'minint_cod')

assignMacroRegion <- function(x) {
  if (x %in% c('Piemonte','Lombardia',
               'Liguria',"Valle d'Aosta/Vallée d'Aoste")) {
    return('North-west')
  }
  if(x %in% c('Veneto','Friuli-Venezia Giulia',
              "Trentino-Alto Adige/Südtirol",
              "Emilia-Romagna")) {
    return('North-east')
  }
  if (x %in% c("Toscana", "Umbria", "Marche", "Lazio")) {
    return('Centre')
  }
  if (x %in% c("Abruzzo", "Puglia", "Basilicata", 
               "Campania", "Calabria", "Molise")) {
    return('South')
  }
  if (x %in% c("Sicilia", "Sardegna")) {
    return("Islands")
  }
}

res_geo_details <- 
  dbGetQuery(conn, "SELECT minint_geopol_cod AS minint_cod,
             ST_Area(ST_Transform(geom, 32632)) AS area
             FROM ita_ge18.subcomuni2018_camera_geo;")

res_demographics <-
  dbGetQuery(conn, "SELECT * FROM ita_ge18.subcomuni2018_census2011_demographics;")

camera_sp$macro_region <- 
  sapply(camera_sp$regione, assignMacroRegion)
camera_sp$macro_region <- 
  factor(camera_sp$macro_region, 
         levels = c("North-west", "North-east",
                    "Centre", "South", "Islands" ))

camera_sp <- 
  merge(camera_sp, res_geo_details, by = 'minint_cod')
camera_sp <- 
  merge(camera_sp, res_demographics, by.x = 'minint_cod', by.y = 'minint_geopol_cod')

camera_sp$pop_density <- 
  camera_sp$pop2011/camera_sp$area
camera_sp$unemployment <- 
  camera_sp$popunemployed2011/camera_sp$popworkforce2011

camera_sp$partwf_perc <- 
  camera_sp$popworkforce2011 / 
  (camera_sp$popworkforce2011 + camera_sp$popnotworkforce2011)

camera_sp$housewive_perc <- 
  camera_sp$pophousewife2011/camera_sp$popworkforce2011
camera_sp$over65_perc <- 
  camera_sp$pop65over2011/camera_sp$pop2011
camera_sp$degree_perc <- 
  camera_sp$popdegree2011/camera_sp$pop2011

camera_sp$foreignpop_perc <- 
  camera_sp$foreignpop2011/camera_sp$pop2011
camera_sp$foreignpop_africa_perc <- 
  camera_sp$foreignpopafrica2011/camera_sp$pop2011

camera_sp$lon <- coordinates(camera_sp)[,1]
camera_sp$lat <- coordinates(camera_sp)[,2]


# Meetup 

## Density

ita_sp <- 
  readOGR('/Users/francesco/Desktop/GIS_Data/Administrative units/Italy/ITA_adm/',
          'ITA_adm0_simp')
ita_rast <- raster(ita_sp)
res(ita_rast) <- 0.2
ita_rast <- rasterize(ita_sp, ita_rast)
quads <- as(ita_rast, 'SpatialPolygons')

# Mean size of square -> mean(sqrt(area(ita_rast)@data@values))
# 19.22603 km

## 2018
meetup2018_event <- 
  dbReadTable(conn, c('ita_ge18','meetup2018_event'))
meetup2018_event_all <- 
  meetup2018_event
meetup2018_comuni2018 <- 
  dbReadTable(conn, c('ita_ge18','meetup2018_comuni2018'))
meetup2018_venue <- 
  dbGetQuery(conn, "SELECT venue_id, lon, lat FROM ita_ge18.meetup2018_venue")

meetup2018_event <- 
  merge(meetup2018_event, meetup2018_comuni2018, by.x = 'venue', by.y='venue_id', all = FALSE)

# Double counting 
meetup2018_event <- 
  meetup2018_event[!duplicated(meetup2018_event[,c("venue","time")]),]

# Remove zero rsvp
meetup2018_event <- 
  meetup2018_event[meetup2018_event$yes_rsvp_count>0,]

meetup2018_event$time <- 
  as.Date(meetup2018_event$time)
meetup2018_event <- 
  subset(meetup2018_event, 
         time >= as.Date('2018-03-04')-90 &
           time <= as.Date('2018-03-04'))

meetup2018_event <- 
  merge(meetup2018_event, meetup2018_venue, by.x = 'venue', by.y = 'venue_id')

meetup2018_event <- 
  SpatialPointsDataFrame(meetup2018_event[c('lon','lat')], 
                         data = meetup2018_event,
                         proj4string = CRS('+proj=longlat +datum=WGS84 +no_defs'))

## 2014
meetup2014_event <- 
  dbReadTable(conn, c('ita_ge18','meetup2014_event'))
meetup2014_event_all <- 
  meetup2014_event
meetup2014_comuni2018 <- 
  dbReadTable(conn, c('ita_ge18','meetup2014_comuni2018'))
meetup2014_venue <- 
  dbGetQuery(conn, "SELECT venue_id, lon, lat FROM ita_ge18.meetup2014_venue")

meetup2014_event <- 
  merge(meetup2014_event, meetup2014_comuni2018, by.x = 'venue', by.y='venue_id', all = FALSE)

# Double counting 
meetup2014_event <- 
  meetup2014_event[!duplicated(meetup2014_event[,c("venue","time")]),]

# Remove zero rsvp
meetup2014_event <- 
  meetup2014_event[meetup2014_event$yes_rsvp_count>0,]

meetup2014_event$time <- 
  as.Date(meetup2014_event$time)
meetup2014_event <- 
  subset(meetup2014_event, 
                           time >= as.Date('2013-02-24')-90 &
                             time <= as.Date('2013-02-24'))
meetup2014_event <- 
  merge(meetup2014_event, meetup2014_venue, 
        by.x = 'venue', by.y = 'venue_id')

meetup2014_event <- 
  SpatialPointsDataFrame(meetup2014_event[
    c('lon','lat')], data = meetup2014_event,
    proj4string = CRS('+proj=longlat +datum=WGS84 +no_defs'))

## Number of recurrent events
utm_crs <- '+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs'
lonlat_crs <- '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'

### 2018
meetup2018_event_utm <- 
  spTransform(meetup2018_event, utm_crs)
meetup2018_event_utm <- 
  meetup2018_event_utm[meetup2018_event_utm$yes_rsvp_count>0,]
meetup2018_event_utm <- 
  meetup2018_event_utm[meetup2018_event_utm$time > as.Date('2018-01-06'),]
meetup2018_event_utm <- 
  meetup2018_event_utm[order(meetup2018_event_utm$time),]

week_seq <- seq(from = as.Date('2018-01-06'),
                to = as.Date('2018-03-04'),
                by = 7)

month_seq <- seq(from = as.Date('2018-01-06'),
                 to = as.Date('2018-03-04'),
                 by = 30)

### Camera
camera_sp_utm <- spTransform(camera_sp, utm_crs)

camera_sp$meetup2018_0km_buff_freq_7days <- 0
camera_sp$meetup2018_0km_buff_freq_30days <- 0

camera_sp$meetup2018_5km_buff_freq_7days <- 0
camera_sp$meetup2018_5km_buff_freq_30days <- 0

camera_sp$meetup2018_10km_buff_freq_7days <- 0
camera_sp$meetup2018_10km_buff_freq_30days <- 0

camera_sp$meetup2018_19km_buff_freq_7days <- 0
camera_sp$meetup2018_19km_buff_freq_30days <- 0

for (i in 1:nrow(camera_sp)) {
  print(i)
  this_poly <- camera_sp_utm[i,]
  this_poly_buff_5k <- gBuffer(this_poly, width = 5000)
  this_poly_buff_10k <- gBuffer(this_poly, width = 10000)
  this_poly_buff_19k <- gBuffer(this_poly, width = 19226)
  
  these_events_0 <- meetup2018_event_utm[this_poly,]
  these_events_5 <- meetup2018_event_utm[this_poly_buff_5k,]
  these_events_10 <- meetup2018_event_utm[this_poly_buff_10k,]
  these_events_19 <- meetup2018_event_utm[this_poly_buff_19k,]
  
  if (nrow(these_events_0) > 0) {
    camera_sp$meetup2018_0km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_0$time, breaks = week_seq)))
    camera_sp$meetup2018_0km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_0$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_5) > 0) {
    camera_sp$meetup2018_5km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_5$time, breaks = week_seq)))
    camera_sp$meetup2018_5km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_5$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_10) > 0) {
    camera_sp$meetup2018_10km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_10$time, breaks = week_seq)))
    camera_sp$meetup2018_10km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_10$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_19) > 0) {
    camera_sp$meetup2018_19km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_19$time, breaks = week_seq)))
    camera_sp$meetup2018_19km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_19$time, breaks = month_seq)))
  }
  
}

### Comuni
comuni_sp_utm <- spTransform(comuni_sp, utm_crs)

comuni_sp$meetup2018_0km_buff_freq_7days <- 0
comuni_sp$meetup2018_0km_buff_freq_30days <- 0

comuni_sp$meetup2018_5km_buff_freq_7days <- 0
comuni_sp$meetup2018_5km_buff_freq_30days <- 0

comuni_sp$meetup2018_10km_buff_freq_7days <- 0
comuni_sp$meetup2018_10km_buff_freq_30days <- 0

comuni_sp$meetup2018_19km_buff_freq_7days <- 0
comuni_sp$meetup2018_19km_buff_freq_30days <- 0

for (i in 1:nrow(comuni_sp)) {
  print(i)
  this_poly <- comuni_sp_utm[i,]
  this_poly_buff_5k <- gBuffer(this_poly, width = 5000)
  this_poly_buff_10k <- gBuffer(this_poly, width = 10000)
  this_poly_buff_19k <- gBuffer(this_poly, width = 19226)
  
  these_events_0 <- meetup2018_event_utm[this_poly,]
  these_events_5 <- meetup2018_event_utm[this_poly_buff_5k,]
  these_events_10 <- meetup2018_event_utm[this_poly_buff_10k,]
  these_events_19 <- meetup2018_event_utm[this_poly_buff_19k,]
  
  if (nrow(these_events_0) > 0) {
    comuni_sp$meetup2018_0km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_0$time, breaks = week_seq)))
    comuni_sp$meetup2018_0km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_0$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_5) > 0) {
    comuni_sp$meetup2018_5km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_5$time, breaks = week_seq)))
    comuni_sp$meetup2018_5km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_5$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_10) > 0) {
    comuni_sp$meetup2018_10km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_10$time, breaks = week_seq)))
    comuni_sp$meetup2018_10km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_10$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_19) > 0) {
    comuni_sp$meetup2018_19km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_19$time, breaks = week_seq)))
    comuni_sp$meetup2018_19km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_19$time, breaks = month_seq)))
  }
  
}

### 2014
meetup2014_event_utm <- 
  spTransform(meetup2014_event, utm_crs)
meetup2014_event_utm <- 
  meetup2014_event_utm[meetup2014_event_utm$yes_rsvp_count>0,]
meetup2014_event_utm <- 
  meetup2014_event_utm[meetup2014_event_utm$time > as.Date('2013-01-06'),]
meetup2014_event_utm <- 
  meetup2014_event_utm[order(meetup2014_event_utm$time),]

week_seq <- seq(from = as.Date('2013-01-06'),
                to = as.Date('2013-02-24'),
                by = 7)

month_seq <- seq(from = as.Date('2013-01-06'),
                 to = as.Date('2013-02-24'),
                 by = 30)

### Camera
camera_sp_utm <- spTransform(camera_sp, utm_crs)

camera_sp$meetup2014_0km_buff_freq_7days <- 0
camera_sp$meetup2014_0km_buff_freq_30days <- 0

camera_sp$meetup2014_5km_buff_freq_7days <- 0
camera_sp$meetup2014_5km_buff_freq_30days <- 0

camera_sp$meetup2014_10km_buff_freq_7days <- 0
camera_sp$meetup2014_10km_buff_freq_30days <- 0

camera_sp$meetup2014_19km_buff_freq_7days <- 0
camera_sp$meetup2014_19km_buff_freq_30days <- 0

for (i in 1:nrow(camera_sp)) {
  print(i)
  this_poly <- camera_sp_utm[i,]
  this_poly_buff_5k <- gBuffer(this_poly, width = 5000)
  this_poly_buff_10k <- gBuffer(this_poly, width = 10000)
  this_poly_buff_19k <- gBuffer(this_poly, width = 19226)
  
  these_events_0 <- meetup2014_event_utm[this_poly,]
  these_events_5 <- meetup2014_event_utm[this_poly_buff_5k,]
  these_events_10 <- meetup2014_event_utm[this_poly_buff_10k,]
  these_events_19 <- meetup2014_event_utm[this_poly_buff_19k,]
  
  if (nrow(these_events_0) > 0) {
    camera_sp$meetup2014_0km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_0$time, breaks = week_seq)))
    camera_sp$meetup2014_0km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_0$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_5) > 0) {
    camera_sp$meetup2014_5km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_5$time, breaks = week_seq)))
    camera_sp$meetup2014_5km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_5$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_10) > 0) {
    camera_sp$meetup2014_10km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_10$time, breaks = week_seq)))
    camera_sp$meetup2014_10km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_10$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_19) > 0) {
    camera_sp$meetup2014_19km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_19$time, breaks = week_seq)))
    camera_sp$meetup2014_19km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_19$time, breaks = month_seq)))
  }
  
}

### Comuni
comuni_sp_utm <- spTransform(comuni_sp, utm_crs)

comuni_sp$meetup2014_0km_buff_freq_7days <- 0
comuni_sp$meetup2014_0km_buff_freq_30days <- 0

comuni_sp$meetup2014_5km_buff_freq_7days <- 0
comuni_sp$meetup2014_5km_buff_freq_30days <- 0

comuni_sp$meetup2014_10km_buff_freq_7days <- 0
comuni_sp$meetup2014_10km_buff_freq_30days <- 0

comuni_sp$meetup2014_19km_buff_freq_7days <- 0
comuni_sp$meetup2014_19km_buff_freq_30days <- 0

for (i in 1:nrow(comuni_sp)) {
  print(i)
  this_poly <- comuni_sp_utm[i,]
  this_poly_buff_5k <- gBuffer(this_poly, width = 5000)
  this_poly_buff_10k <- gBuffer(this_poly, width = 10000)
  this_poly_buff_19k <- gBuffer(this_poly, width = 19226)
  
  these_events_0 <- meetup2014_event_utm[this_poly,]
  these_events_5 <- meetup2014_event_utm[this_poly_buff_5k,]
  these_events_10 <- meetup2014_event_utm[this_poly_buff_10k,]
  these_events_19 <- meetup2014_event_utm[this_poly_buff_19k,]
  
  if (nrow(these_events_0) > 0) {
    comuni_sp$meetup2014_0km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_0$time, breaks = week_seq)))
    comuni_sp$meetup2014_0km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_0$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_5) > 0) {
    comuni_sp$meetup2014_5km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_5$time, breaks = week_seq)))
    comuni_sp$meetup2014_5km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_5$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_10) > 0) {
    comuni_sp$meetup2014_10km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_10$time, breaks = week_seq)))
    comuni_sp$meetup2014_10km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_10$time, breaks = month_seq)))
  }
  
  if (nrow(these_events_19) > 0) {
    comuni_sp$meetup2014_19km_buff_freq_7days[i] <- 
      length(unique(cut(these_events_19$time, breaks = week_seq)))
    comuni_sp$meetup2014_19km_buff_freq_30days[i] <- 
      length(unique(cut(these_events_19$time, breaks = month_seq)))
  }
  
}


## Actual density

### 2018
meetup2018_event_r <- 
  rasterize(coordinates(meetup2018_event), ita_rast, fun='count', background=0)
meetup2018_event_spdf <- as(meetup2018_event_r, "SpatialPixelsDataFrame")
meetup2018_event_df <- as.data.frame(meetup2018_event_spdf)
colnames(meetup2018_event_df) <- c("value", "x", "y")

v2018 <- extract(meetup2018_event_r, comuni_sp, fun = mean)
comuni_meetup2018_density <- data.frame(v2018)
comuni_meetup2018_density$v2018[is.na(comuni_meetup2018_density$v2018)] <- 0
comuni_sp$meetup2018_19km_quad_density <- comuni_meetup2018_density$v2018

v2018 <- extract(meetup2018_event_r, camera_sp, fun = mean)
camera_meetup2018_density <- data.frame(v2018)
camera_meetup2018_density$v2018[is.na(camera_meetup2018_density$v2018)] <- 0
camera_sp$meetup2018_19km_quad_density <- camera_meetup2018_density$v2018

### 2014
meetup2014_event_r <- 
  rasterize(coordinates(meetup2014_event), ita_rast, fun='count', background=0)
meetup2014_event_spdf <- as(meetup2014_event_r, "SpatialPixelsDataFrame")
meetup2014_event_df <- as.data.frame(meetup2014_event_spdf)
colnames(meetup2014_event_df) <- c("value", "x", "y")

v2014 <- extract(meetup2014_event_r, comuni_sp, fun = mean)
comuni_meetup2014_density <- data.frame(v2014)
comuni_meetup2014_density$v2014[is.na(comuni_meetup2014_density$v2014)] <- 0
comuni_sp$meetup2014_19km_quad_density <- comuni_meetup2014_density$v2014

v2014 <- extract(meetup2014_event_r, camera_sp, fun = mean)
camera_meetup2014_density <- data.frame(v2014)
camera_meetup2014_density$v2014[is.na(camera_meetup2014_density$v2014)] <- 0
camera_sp$meetup2014_19km_quad_density <- camera_meetup2014_density$v2014

## Density (Hexagon)

library(raster)  
library(rgdal)
ita_sp <- 
  readOGR('/Users/francesco/Desktop/GIS_Data/Administrative units/Italy/ITA_adm/',
          'ITA_adm0_simp')
ita_sp <- spTransform(ita_sp, CRS("+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"))

ita_sp_buff_10km <- rgeos::gBuffer(ita_sp, width = 10000)
ita_sp_buff_30km <- rgeos::gBuffer(ita_sp, width = 30000)

size <- 10000
hex_points_10km <- spsample(ita_sp_buff_10km, type = "hexagonal", cellsize = size)
hex_grid_10km <- HexPoints2SpatialPolygons(hex_points_10km, dx = size)

size <- 25000
hex_points_25km <- spsample(ita_sp_buff_10km, type = "hexagonal", cellsize = size)
hex_grid_25km <- HexPoints2SpatialPolygons(hex_points_25km, dx = size)

size <- 50000
hex_points_50km <- spsample(ita_sp_buff_30km, type = "hexagonal", cellsize = size)
hex_grid_50km <- HexPoints2SpatialPolygons(hex_points_50km, dx = size)

require(rgeos)
hex_grid_10km <- 
  gIntersection(hex_grid_10km, ita_sp, byid = TRUE)
hex_grid_25km <- 
  gIntersection(hex_grid_25km, ita_sp, byid = TRUE)
hex_grid_50km <- 
  gIntersection(hex_grid_50km, ita_sp, byid = TRUE)

hex_grid_10km <- 
  spTransform(hex_grid_10km, CRS("+proj=longlat +datum=WGS84 +no_defs"))
hex_grid_25km <- 
  spTransform(hex_grid_25km, CRS("+proj=longlat +datum=WGS84 +no_defs"))
hex_grid_50km <- 
  spTransform(hex_grid_50km, CRS("+proj=longlat +datum=WGS84 +no_defs"))

proj4string(meetup2018_event) <- 
  CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj4string(meetup2014_event) <- 
  CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

hex_grid_10km <- 
  SpatialPolygonsDataFrame(hex_grid_10km, 
                           data.frame(id = 1:length(hex_grid_10km)),
                           match.ID = FALSE)

hex_grid_25km <- 
  SpatialPolygonsDataFrame(hex_grid_25km, 
                           data.frame(id = 1:length(hex_grid_25km)),
                           match.ID = FALSE)

hex_grid_50km <- 
  SpatialPolygonsDataFrame(hex_grid_50km, 
                           data.frame(id = 1:length(hex_grid_50km)),
                           match.ID = FALSE)

hex_grid_10km_density <- 
  as.data.frame(table(over(meetup2018_event, hex_grid_10km)$id))
hex_grid_25km_density <- 
  as.data.frame(table(over(meetup2018_event, hex_grid_25km)$id))
hex_grid_50km_density <- 
  as.data.frame(table(over(meetup2018_event, hex_grid_50km)$id))

hex_grid_10km_meetup2018 <- 
  merge(hex_grid_10km, hex_grid_10km_density, by.x = 'id', by.y = 'Var1')
hex_grid_10km_meetup2018$Freq[
  is.na(hex_grid_10km$Freq_meetup2018)] <- 0

hex_grid_25km_meetup2018 <- 
  merge(hex_grid_25km, hex_grid_25km_density, by.x = 'id', by.y = 'Var1')
hex_grid_25km_meetup2018$Freq[
  is.na(hex_grid_25km_meetup2018$Freq)] <- 0

hex_grid_50km_meetup2018 <- 
  merge(hex_grid_50km, hex_grid_50km_density, by.x = 'id', by.y = 'Var1')
hex_grid_50km_meetup2018$Freq[is.na(hex_grid_50km_meetup2018$Freq)] <- 0

hex_grid_10km_density <- 
  as.data.frame(table(over(meetup2014_event, hex_grid_10km)$id))
hex_grid_25km_density <- 
  as.data.frame(table(over(meetup2014_event, hex_grid_25km)$id))
hex_grid_50km_density <- 
  as.data.frame(table(over(meetup2014_event, hex_grid_50km)$id))

hex_grid_10km_meetup2014 <- merge(hex_grid_10km, hex_grid_10km_density, by.x = 'id', by.y = 'Var1')
hex_grid_10km_meetup2014$Freq[is.na(hex_grid_10km_meetup2014$Freq)] <- 0

hex_grid_25km_meetup2014 <- merge(hex_grid_25km, hex_grid_25km_density, by.x = 'id', by.y = 'Var1')
hex_grid_25km_meetup2014$Freq[is.na(hex_grid_25km_meetup2014$Freq)] <- 0

hex_grid_50km_meetup2014 <- merge(hex_grid_50km, hex_grid_50km_density, by.x = 'id', by.y = 'Var1')
hex_grid_50km_meetup2014$Freq[is.na(hex_grid_50km_meetup2014$Freq)] <- 0


## Load simplified geom and add data to postgis

comuni_sp_keep005 <- 
  readOGR("scrape_2018_ita_ge_results/output_shapefiles/",
          'istat_comuni_2018_005')
comuni_sp_keep005 <- st_as_sf(comuni_sp_keep005)
comuni_to_merge <- c("001130", "001265", "001063", "001028", "001316")
comuni_sp_keep005$PRO_COM_T <- as.character(comuni_sp_keep005$PRO_COM_T)
comuni_sp_keep005$PRO_COM_T[comuni_sp_keep005$PRO_COM_T %in% comuni_to_merge] <- "000000"
comuni_sp_keep005 <- comuni_sp_keep005[,c('PRO_COM_T')]
comuni_sp_keep005 <-
  comuni_sp_keep005 %>%
  dplyr::group_by(PRO_COM_T) %>%
  dplyr::summarise(STATUS = 'dissolved')
comuni_sp_keep005$STATUS <- NULL
comuni_sp_keep005 <- as(comuni_sp_keep005, 'Spatial')

camera_sp_keep005 <- 
  readOGR("scrape_2018_ita_ge_results/output_shapefiles/",
          'minint_camera_005')

# # Neighborhood
library(spdep)
camera_nb <- list()
comuni_nb <- list()

for(this_km in c(10, 25, 50)) {
  camera_nb[[as.character(this_km)]] <-
    dnearneigh(coordinates(camera_sp),  
               d1=0, d2=this_km, longlat=T, row.names=camera_sp$minint_cod)
  comuni_nb[[as.character(this_km)]] <-
    dnearneigh(coordinates(comuni_sp),  
               d1=0, d2=this_km, longlat=T, row.names=comuni_sp$PRO_COM_T)
}

# Remoteness
IT <- read.delim("~/Downloads/IT (2)/IT.txt", header=FALSE)
IT <- subset(IT, grepl("^PPL", V8) & V15 > 9999)

IT_point <- 
  SpatialPointsDataFrame(coords = IT[,c('V6','V5')],
                         data = IT, 
                         proj4string = 
                           CRS('+proj=longlat +datum=WGS84 +no_defs'))

utmStr <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
comuni_sp_utm <- spTransform(comuni_sp, utmStr)
camera_sp_utm <- spTransform(camera_sp, utmStr)
IT_point <- spTransform(IT_point, utmStr)

comuni_sp$dist_10k_ppl <- NA
for (i in 1:nrow(comuni_sp)) {
  print(i)
  res <- gDistance(comuni_sp_utm[i,], IT_point, byid = TRUE)
  comuni_sp$dist_10k_ppl[i] <-  min(res[,1])
}
camera_sp$dist_10k_ppl <- NA
for (i in 1:nrow(camera_sp)) {
  print(i)
  res <- gDistance(camera_sp_utm[i,], IT_point, byid = TRUE)
  camera_sp$dist_10k_ppl[i] <-  min(res[,1])
}

save(camera_sp, comuni_sp, 
     camera_sp_keep005, comuni_sp_keep005,
     camera_nb, comuni_nb,
     meetup2014_event, meetup2018_event,
     hex_grid_10km_meetup2018, hex_grid_25km_meetup2018, hex_grid_50km_meetup2018,
     hex_grid_10km_meetup2014, hex_grid_25km_meetup2014, hex_grid_50km_meetup2014,
     file = '00_prepare_data.RData')

